/**
 * @file Mesh.h
 * @author Lishijie (lsj1018845759@outlook.com)
 * @brief 
 * @version 0.1
 * @date 2020-11-08
 * 
 * @copyright Copyright (c) 2020
 * 
 */
#ifndef _MESH_H_
#define _MESH_H_
#include "Domain.h"
#include "Dofs.h"

#define TEMPLATE template<Dimension DIM>
/**
 * @brief 有限元计算网格类，所有网格从这个类派生；
 *        功能：设置和返回网格大小，计算并返回指定全局编号对应的自由度信息，
 *             返回指定单元的所有自由度的全局编号，以及其他网格信息；
 */
TEMPLATE 
class Mesh
{
public:
    /**
     * @brief 该网格中单元个数；
     * 
     * @return long int 
     */
    virtual long int n_element() = 0;
    /**
     * @brief 该网格中节点个数;
     * 
     * @return long int 
     */
    virtual long int n_point() = 0; 
    /**
     * @brief x坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& x_seg() = 0;
    /**
     * @brief y坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& y_seg() = 0;
    /**
     * @brief z坐标轴上划分的段数；
     * 
     * @return long& 
     */
    virtual long int& z_seg() = 0;
    /**
     * @brief x坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real x_h() = 0;
    /**
     * @brief y坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real y_h() = 0;
    /**
     * @brief z坐标轴上每一段长度；
     * 
     * @return long& 
     */
    virtual Real z_h() = 0;
    /**
     * @brief 返回全局编号为i的自由度
     * 
     * @param i 指定全局编号
     * @return Dofs<DIM> 对应的自由度；
     */
    virtual Dofs<DIM> DofsofIndex(int i) = 0;
    /**
     * @brief 返回指定单元所有全局编号；
     * 
     * @param i 指定单元序数；
     * @return std::vector<int> 存放该单元所有全局编号的数组
     */
    virtual std::vector<int> NodeofEle(int i) = 0;
    /**
     * @brief 返回所有边界上点的全局编号；
     * 
     * @return std::vector<int>& 存放该网格所有边界上自由度全局编号的数组
     */
    virtual std::vector<int>& Boundary() = 0;
    /**
     * @brief 返回计算区域的指针
     * 
     * @return RectangleDomain* 计算区域的指针
     */
    virtual RectangleDomain* RDomain() = 0;
    /**
     * @brief 二维度情况下，第i,j号自由度对应的全局编号
     * 
     * @param i 下标１
     * @param j 下标２
     * @return int 全局编号
     */
    virtual int CortoIdx(int i, int j) = 0;
    /**
     * @brief 判断该index序号点是否一个边界点； 
     * 
     * @param index 需要判断点的全局编号
     * @return true 是边界点
     * @return false 不是边界点
     */
    virtual bool IsBoundary(int index) = 0;
protected:
};
#undef TEMPLATE
/**
 * @brief Q1网格，对应单元为Q1单元；
 * 
 */
/**********************************************************/
/*********************** Q1 2D ****************************/
/**********************************************************/

class Q1Mesh:public Mesh<2>
{
public:
    Q1Mesh(){};
    Q1Mesh(RectangleDomain*);
    Q1Mesh(RectangleDomain*,long int[]);
    Q1Mesh(RectangleDomain*,std::vector<long int>);
    Q1Mesh(std::shared_ptr<RectangleDomain>);
    Q1Mesh(std::shared_ptr<RectangleDomain>,long int[]);
    Q1Mesh(std::shared_ptr<RectangleDomain>,std::vector<long int>);
    ~Q1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_point(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return abs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) /static_cast<Real>(x_seg());};
    virtual Real y_h(){return abs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) /static_cast<Real>(y_seg());};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};//二维情况不用，接口就私有吧
    virtual Real z_h(){return y_h();};//二维情况用不到，接口就私有吧
    void ReBoundary();
};
Q1Mesh::Q1Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

Q1Mesh::Q1Mesh(RectangleDomain* Domain,long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg,seg+2,_seg);
    ReBoundary();
}

Q1Mesh::Q1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
}

Q1Mesh::Q1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
}

long int Q1Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q1Mesh::n_point()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1));
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q1Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (x_seg() + 1);
    int j = index / static_cast<int>(x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h();
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q1Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (x_seg());
    int b = i % (x_seg());
    Idx.push_back(a * (x_seg() + 1) + b);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (x_seg() + 1));
    Idx.push_back(Idx[2] - 1);
    return Idx;
}
void Q1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_point()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
}
/**
 * @brief P1网格，对应单元为P1单元；
 * 
 */
/**********************************************************/
/*********************** P1 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第2个顶点连线，划分为两个三角形，从左至右排序。
 * 2  ——  3 
 * |     /|    如图左所示，0-3-2, 0-1-3 构成两个 P1 单元片，
 * |    / |    且 0-3-2 编号在前。
 * |   /  |
 * |  /   |
 * | /    |
 * |/     |
 * 0  ——  1
 *
 */

class P1Mesh:public Mesh<2>
{
public:
    P1Mesh(){};
    P1Mesh(RectangleDomain*);
    P1Mesh(RectangleDomain*, long int[]);
    P1Mesh(RectangleDomain*, std::vector<long int>);
    P1Mesh(std::shared_ptr<RectangleDomain>);
    P1Mesh(std::shared_ptr<RectangleDomain>, long int[]);
    P1Mesh(std::shared_ptr<RectangleDomain>, std::vector<long int>);
    ~P1Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_point(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return abs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) / static_cast<Real>(x_seg());};
    virtual Real y_h(){return abs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) / static_cast<Real>(y_seg());};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
        virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};
    virtual Real z_h(){return y_h();};
    void ReBoundary();
};

P1Mesh::P1Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

P1Mesh::P1Mesh(RectangleDomain* Domain, long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg, seg+2, _seg);
    ReBoundary();
}

P1Mesh::P1Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
}

P1Mesh::P1Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
}

long int P1Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}

long int P1Mesh::n_point()
{
    return static_cast<long int>((_seg[0] + 1) * (_seg[1] + 1)) ;
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> P1Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (x_seg() + 1);
    int j = index / static_cast<int>(x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h();
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h();
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> P1Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (2 * x_seg());
    int b = i % (2 * x_seg());
    if (b % 2 == 0)
    {
	Idx.push_back(a * (x_seg() + 1) + b / 2);
	Idx.push_back(Idx[0] + x_seg() + 2);
	Idx.push_back(Idx[1] - 1);
    }
    else
    {
	Idx.push_back(a * (x_seg() + 1) + b / 2);
	Idx.push_back(Idx[0] + 1);
	Idx.push_back(Idx[1] + x_seg() + 1);
    }
    return Idx;
}

void P1Mesh::ReBoundary()
{
    for(int i = 0;i <= _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_point()- i - 1);
    }
    for(int j = 1;j <= _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j*(_seg[0] + 1));
        _BndIndex.push_back((j + 1)*(_seg[0] + 1) - 1);
    }
}
/**
 * @brief Q2网格，对应单元为Q2单元；
 * 
 */
/**********************************************************/
/*********************** P2 2D ****************************/
/**********************************************************/
class Q2Mesh:public Mesh<2>
{
public:
    Q2Mesh(){};
    Q2Mesh(RectangleDomain*);
    Q2Mesh(RectangleDomain*,long int[]);
    Q2Mesh(RectangleDomain*,std::vector<long int>);
    Q2Mesh(std::shared_ptr<RectangleDomain>);
    Q2Mesh(std::shared_ptr<RectangleDomain>,long int[]);
    Q2Mesh(std::shared_ptr<RectangleDomain>,std::vector<long int>);
    ~Q2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_point(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return abs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) /static_cast<Real>(x_seg());};
    virtual Real y_h(){return abs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) /static_cast<Real>(y_seg());};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};//二维情况不用，接口就私有吧
    virtual Real z_h(){return y_h();};//二维情况用不到，接口就私有吧
    void ReBoundary();
};
Q2Mesh::Q2Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

Q2Mesh::Q2Mesh(RectangleDomain* Domain,long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg,seg+2,_seg);
    ReBoundary();
}

Q2Mesh::Q2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
}

Q2Mesh::Q2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
}

long int Q2Mesh::n_element()
{
    return static_cast<long int>(_seg[0] * _seg[1]);
}

long int Q2Mesh::n_point()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1));
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> Q2Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (2 * x_seg() + 1);
    int j = index / static_cast<int>(2 * x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h() /2.0;
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h() /2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
} 

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> Q2Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (x_seg());
    int b = i % (x_seg());
    Idx.push_back(a * 2 * (2 * x_seg() + 1) + 2 * b);
    Idx.push_back(Idx[0] + 2);
    Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
    Idx.push_back(Idx[2] - 2);
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + (2 * x_seg() + 1));
    Idx.push_back(Idx[2] - 1);
    Idx.push_back(Idx[0] + (2 * x_seg() + 1));
    Idx.push_back(Idx[7] + 1);
    return Idx;
}
void Q2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_point()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}
/**
 * @brief P2网格，对应单元为P2单元；
 * 
 */
/**********************************************************/
/*********************** P2 2D ****************************/
/**********************************************************/

/* 如何划分单元格？                     
 * 首先将实际计算区域划分为规则的矩形网格，对于每个矩形网格，
 * 将其第0个和第8个顶点连线，划分为两个三角形，从左至右排序。
 * 从左下角为第一个号，逆时针排序．
 * 6  —7—  8 
 * |     /|    如图左所示，0-8-6-4-7-3, 0-2-8-1-5-4 构成两个 P2 单元片，
 * |    / |    且 0-8-6-4-7-3 编号在前。
 * 3   4  5
 * |  /   |
 * | /    |
 * |/     |
 * 0  —1—  2
 *
 */

class P2Mesh:public Mesh<2>
{
public:
    P2Mesh(){};
    P2Mesh(RectangleDomain*);
    P2Mesh(RectangleDomain*, long int[]);
    P2Mesh(RectangleDomain*, std::vector<long int>);
    P2Mesh(std::shared_ptr<RectangleDomain>);
    P2Mesh(std::shared_ptr<RectangleDomain>, long int[]);
    P2Mesh(std::shared_ptr<RectangleDomain>, std::vector<long int>);
    ~P2Mesh(){delete _RDomain;};
    virtual long int n_element();
    virtual long int n_point(); 
    virtual long int& x_seg(){return _seg[0];};
    virtual long int& y_seg(){return _seg[1];}; 
    virtual Real x_h(){return abs(_RDomain->LeftBottomPoint()[0] - _RDomain->RightBottomPoint()[0]) / static_cast<Real>(x_seg());};
    virtual Real y_h(){return abs(_RDomain->LeftBottomPoint()[1] - _RDomain->LeftTopPoint()[1]) / static_cast<Real>(y_seg());};
    virtual Dofs<2> DofsofIndex(int i);
    virtual std::vector<int> NodeofEle(int i);
    virtual std::vector<int>& Boundary(){return _BndIndex;};
    virtual RectangleDomain* RDomain() {return _RDomain;};
    virtual int CortoIdx(int i, int j) {return (i*(_seg[0]+1) + j);};
    virtual bool IsBoundary(int index){return !(std::find(_BndIndex.begin(),_BndIndex.end(),index)==_BndIndex.end());};
private:
    RectangleDomain* _RDomain;
    long int _seg[2];
    std::vector<int> _BndIndex;
    virtual long int& z_seg(){return _seg[1];};
    virtual Real z_h(){return y_h();};
    void ReBoundary();
};

P2Mesh::P2Mesh(RectangleDomain* Domain)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

P2Mesh::P2Mesh(RectangleDomain* Domain, long int seg[])
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg, seg+2, _seg);
    ReBoundary();
}

P2Mesh::P2Mesh(RectangleDomain* Domain,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(Domain);
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();

}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    _seg[0] = 5;
    _seg[1] = 5;
    ReBoundary();
}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,long int seg[])
{
     _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
     std::copy(seg,seg+2,_seg);
     ReBoundary();
}

P2Mesh::P2Mesh(std::shared_ptr<RectangleDomain> _Domain_sptr,std::vector<long int> seg)
{
    _RDomain = static_cast<RectangleDomain*>(_Domain_sptr.get());
    std::copy(seg.begin(),seg.end(),_seg);
    ReBoundary();
}

long int P2Mesh::n_element()
{
    return static_cast<long int>(2 * _seg[0] * _seg[1]);
}

long int P2Mesh::n_point()
{
    return static_cast<long int>((2 * _seg[0] + 1) * (2 * _seg[1] + 1)) ;
}

/** 
 * @brief 返回全局编号对应的自由度
 * 
 * @param index 
 * 
 * @return 
 */
Dofs<2> P2Mesh::DofsofIndex(int index)
{
    Dofs<2> temp_dofs;
    int i = index % (2 * x_seg() + 1);
    int j = index / static_cast<int>(2 * x_seg() + 1);
    Real x_c = _RDomain->LeftBottomPoint()[0] + i * x_h() / 2.0;
    Real y_c = _RDomain->LeftBottomPoint()[1] + j * y_h() / 2.0;
    temp_dofs[0] = x_c;
    temp_dofs[1] = y_c;
    temp_dofs.GlobalIndex() = index;
    return temp_dofs;
}

/** 
 * @brief 按顺序返回指定单元包含的所有自由度的编号
 * 
 * @param i 
 * 
 * @return 
 */
std::vector<int> P2Mesh::NodeofEle(int i)
{
    std::vector<int> Idx;
    int a = i / (2 * x_seg());
    int b = i % (2 * x_seg());
    if (b % 2 == 0)
    {
	Idx.push_back(a * 2 * (2 * x_seg() + 1) + b);
	Idx.push_back(Idx[0] + 2 * (2 * x_seg() + 1) + 2);
	Idx.push_back(Idx[1] - 2);
    Idx.push_back(Idx[0] + 1 + (2 * x_seg() + 1));
    Idx.push_back(Idx[1] - 1);
    Idx.push_back(Idx[3] - 1);
    }
    else
    {
	Idx.push_back(a * 2 *(2 * x_seg() + 1) + b - 1);
	Idx.push_back(Idx[0] + 2);
	Idx.push_back(Idx[1] + 2 * (2 * x_seg() + 1));
    Idx.push_back(Idx[0] + 1);
    Idx.push_back(Idx[1] + 2 *x_seg() + 1);
    Idx.push_back(Idx[4] - 1);
    }
    return Idx;
}

void P2Mesh::ReBoundary()
{
    for(int i = 0;i <= 2 * _seg[0];i++)
    {
        _BndIndex.push_back(i);
        _BndIndex.push_back(n_point()- i - 1);
    }
    for(int j = 1;j <= 2 * _seg[1] - 1;j++)
    {
        _BndIndex.push_back(j* (2  * _seg[0] + 1));
        _BndIndex.push_back((j + 1)*(2 * _seg[0] + 1) - 1);
    }
}
#else
#endif
